Radial Basis Function Network (rbfn) Approximation of Finite Element Models for Real-time Simulation

Nonlinearities inherent in soft-tissue interactions create roadblocks to realization of high-fidelity real-time haptics-based medical simulations. While finite element (FE) formulations offer greater accuracy over conventional spring-mass-network models, computational-complexity limits achievable simulation-update rates. Direct interaction with sensorized physical surrogates, in offline or online modes, allows a temporary sidestepping of computational issues but hinders parametric analysis and true exploitation of a simulation-based testing paradigm. Hence, in this paper, we develop Radial-Basis Neural-Network approximations, to FE-model data within a Modified Resource Allocating Network (MRAN) framework. Real-time simulation of the reduced order neural-network approximations at high temporal resolution provided the haptic-feedback. Validation studies are being conducted to evaluate the kinesthetic realism of these models with medical experts. INTRODUCTION Virtual reality (VR) simulation based training holds immense promise for offline development of skills in arenas ranging from flight simulation to medical procedural training. The promise lies in the ability to not only present controlled sets of stimuli but to simultaneously train both cognitive and sensorimotor skills of the user. One aspect of this work relies on development of haptic user interfaces (HUIs) capable of rendering the forces to the user, which is discussed elsewhere [1]. The other aspect focuses on development of computational haptic models capable of rendering the high fidelity haptic interactions and is discussed here. Over the years VR simulators for medical procedural training have gradually transitioned from a primarily unidirectional visual engagement process to a more bidirectional kinesthetic immersion [2]. Effective implementation of the force-feedback algorithms thus requires high fidelity models which are capable of being computed within the timing constraints prescribed by a deterministic real-time framework necessary for haptics. Enhancing the haptic response (i.e.) to reflect the forceresponse behavior of real bio-physical systems (tissues, organs etc.) makes nonlinear models necessary. Current haptics models of soft tissues (used in COTS medical simulators) rely on relatively simple/ simplistic linear models (spring-mass-damper systems or linear finite element analysis systems) to compute the feedback at high update rates However, the complexity of computing the higher spatial and temporal fidelity response from such nonlinear models creates challenges. While methods such as finite element methods (FEM) offer potential means of computing such responses, they can be pursued only in an offline setting. Hence, in this work, we seek to develop suitable approximation methods to effectively and parametrically capture the physics in reduced order models. The results obtained using FE systems were used as case studies for which radial basis function based approximation model was developed. A modified resource allocating network (MRAN) method was adopted to determine the neural network states and an extended Kalman filter (EKF) method was implemented for optimizing these parameters. The response of the approximation models tend to be faster and easier to simulate than the original nonlinear system. Moreover, methods based on radial basis functions can be trained offline to have an optimal network structure and deployed online in the predictive phase at very high update rates suitable for medical surgical applications. We examine the applicability by studying the approximation of FEM analysis of a linear elastic and plastic cantilever beam. Finally, a virtual-haptic environment was developed using MATLAB Simulink/ VRML to deploy the resulting force reflection models. Working through this simple problem provides a greater insight about the principal issues related to numerical computation of RBFN model as well as implementation of real time haptic interface. BACKGROUND Realistic real-time haptic applications typically require a feedback or update at greater than 30 Hz for visual sensations and minimum of 500 Hz for haptic sensations [3]. Currently VR-and haptic (VR-H) simulations are constrained by both limitations inherent to Haptic-User-Interfaces (HUIs) as well as haptic computational models. However, this is a coupled problem – choices of HUIs determine the sophistication of haptic computational models and the current technology only permits the simplest of models to be run in real-time [4] (i.e., lumped parameter spring/mass systems). While the use of linear FE models helps overcome some of the computational limitations, the accuracy and fidelity of haptic models are always in question. More realistic, complex 3D finite element based soft tissue models, are currently far outside of the realm of real time simulations. Linear elasticity is used for modeling the deformable materials principally due to simplicity of ensuing computations. However, the physical behavior of soft tissue may be considered as linear elastic only for small deformations [5, 6] (typically less than 10% of the mesh size). Thus, linear elastic FE models are not valid for large displacements and are not invariant with respect to rotations [7]. Finally selection of mesh density as well as parameters for the linear elastic systems remains challenging exercise. Thus, currently many surgical haptic simulators depend upon the subject studies to tune these parameters to achieve “realistic performance”. Some of the simpler deformation techniques use surface models, where the masses are concentrated in the mesh vertices connected by springs. Such a surface model contains less mass points than a volumetric model defining the same shape and is therefore more efficient [8]. A surface model, however, is inherently inaccurate and yields physically invalid deformations (e.g. self penetration). Thus, volumetric models are better suited for the simulation of deformable objects, especially when cutting is required. The most commonly used methods for developing soft tissue simulations are the Mass-Spring Method (MSM) and the FEM and method obtained by modification and combination of these two approaches [9]. In all these cases, the choice of the appropriate simulation method is influenced by factors such as computing efficiency, required accuracy and the types of manipulations that have to be performed. With the use of higher fidelity models, the computation time increases rapidly and as a result the real-time haptic performance for surgical simulation has to be compromised upon. In this work, we focus on the issues surrounding numerical implementation of nonlinear models and an alternative novel approach to develop approximation models for the nonlinear FE systems using radial basis function network (RBFN). Such an approach can be made adaptive to different model parameters with nonlinear characteristics while remaining computationally tractable. In addition, using such methods can also prove to be valuable for implementing real-time simulation for other nonlinear FE models outside the medical simulator domain. RADIAL BASIS FUNCTION NETWORKS (RBFN) Over the past few decades, Artificial Neural Networks (ANNs) have emerged as a powerful set of tools in pattern classification, time series analysis, signal processing, dynamical system modeling and control. The popularity of ANNs can be attributed to the fact that these network models are frequently able to learn behavior when traditional modeling is very difficult to generalize. Typically, a neural network consists of several computational nodes called perceptrons arranged in layers. The number of hidden nodes essentially determines the degrees of freedom of the non-parametric model. A small number of hidden units may not be enough to capture a given system’s complex input-output mapping and alternately a large number of hidden units may overfit the data and may not generalize the behavior. An RBF network is a two-layer feed-forward type network in which the input is transformed by the basis functions at the hidden layer. At the output layer, linear combinations of the hidden layer node responses are added to form the output. The name RBF comes from the fact that the basis functions in the hidden layer nodes are radially symmetric. In [11] authors report that the choice of the basis function is not crucial to the performance of the FIG. 1 MULTILAYER HIDDEN NODE NEURAL NETWORK WITH WEIGHTED OUTPUTS [10] network. The most common choice however, is the Gaussian function which can be defined by a mean and a standard deviation. Figure 1 shows a schematic diagram of an RBF network with n, l and m respectively of input, hidden and output layer nodes for the general transformation of nd points of X(x1, x2, x3, ... , xnd) in the input space to points Y(y1, y2, y3, ... , ynd) in the output space. In RBF networks, the connections between the input and the hidden layers are generally not weighted. The inputs therefore reach the hidden layer nodes unchanged. For an input Xi, the j th hidden node produces a response φj given by, where, ||Xi – μj|| is the distance between the point representing the input, Xi and the centre of the j th hidden node, μj as measured by Euclidean norm. R is a state error noise covariance matrix. Generally, it is a positive definite matrix and for circular Gaussian functions, the off diagonal elements of R are zero and diagonal elements are σk 2 , for k = 1...n. The output yik of the network at the output node is given by the weighted sum of RBFs, where, { } 1 2 , , , T T h W w w w − ⋯ is a vector of h linear weights or amplitudes, φ vector of RBFs with centers at μj. In the special case where the number of hidden layer nodes is equal to the number of data in the training set (l = nd) and the RBF centres coincide with the inputs (Uj = Xi, where i = j = 1, 2..., nd), the hidden layer response according to (1) becomes unity for j = i. If the basis functions are truly localized, the response of the other hidden layer nodes will be near zero (i.e. σj for j ≠ i are such that φj ~ 0 for j ≠ i). It can also be seen that (2) gives the exact output when the output layer weight is equal to the output (the contribution to the weighted summation from j = i is yik and that from all j ≠ i is nearly zero). In the ideal case therefore, RBF network can be made to map points in N-dimensional input space exactly on to points in M-dimensional output space. This however, is not practical when nd is large in which case a few input points are chosen to represent the entire input data set. The original RBF method requires that there be as many RBF centres as there are distinct data points in the input space. This however, is not possible in practice because that increases the complexity of the final RBF model tremendously. Moreover, the inputs usually occur in clusters making overlapping of receptive fields inevitable. Choosing all points as RBF centres will therefore lead to more number of redundant nodes with a huge network involving long training and computation times. Modified Resource Allocating Network. The RBFN method that was implemented for our problem is Modified Resource Allocating Network (MRAN). MRAN adopts the basic idea of adaptively “growing” the number of radial basis functions where needed to null local errors, and also includes a “pruning strategy” to eliminate little-needed radial basis functions (those with weights smaller than some tolerance), with the overall goal of finding a minimal RBF network. RAN allocates new units as well as adjusts the network parameters to reflect the complexity of function being approximated. Though the optimal approximation is to add an impulse unit at the data point to match the output error, such a method actually lacks smoothness and is more error prone. Hence, in this work, the method explained in [12] using the Gaussian functions centered at the input data points to achieve the desired output was used. New nodes are not added at every point but are actually restricted by enforcing three main conditions as follows: Equation (3) ensures that a new node is added only if it is sufficiently far from all the existing nodes and equation (4) ensures that only if the approximation error using existing nodes exceeds the error specification. The final condition (equation(5)) compensates for the noise present in the observations data by determining the RMS error of last Nw observations and ensures that a new node is added only if the noise in data exceeds the specified threshold limit. For all the cases to be considered in this work only the standard Gaussian functions are used and the rotational parameters (off diagonal elements of the covariance matrix) are not learned along with the other network parameters, Θ. In this case, if the input space dimension is n, output dimension of the network is 1 and the total number of nodes is h, then the size of the parameter vector considering the individual dimensions of wj, μj and σj for each node, will be: For each observation input to the network, feasibility of the constraints (3), (4) and (5) are determined and a new node will be added if all the conditions hold true. Irrespective of this result, all the network parameters will be updated using Extended Kalman Filter (EKF) as summarized below. Extended Kalman Filter (EKF). This method is used for online adaptation of parameter of our nonlinear function approximation problem. So, the update relationships for each observation input are given below. First, it is necessary to determine the sensitivity (Jacobian) matrix of the RBFN for which a linearized model is used. ( ) 1 2 1 , , , , h T ik j j h j y w W X φ φ μ μ μ = = = ∑ ⋯ (2) ( ) ( ) ( ) 1 exp T j i j i j X R X φ μ μ − = − − −